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INTUOmiCTlON 


In  the  late  lyhO’s  the  aeronomy  branch  at  the  Bit!,  needed  the  solu- 
tions to  sets  of  stiff  ordinary  differential  equations  (01)1  S)  that 
describe  the  jiositive  and  ncjjative  ion  chemistry  in  the  earth's  D-rejtion 
(i0-85  kml  . .\deiiuate  mathematical  techniques  for  handliiij;  stiff  01)1, S 
were  unknown  to  us  at  that  time.  Krejtel,  a jibysicisi,  aiqiroached  this 
problem  empirically  and  developed  a stiff  OIH'.  intejtrator.  This  report 
sketches  the  K-method  of  integration  and  provides  a few  examjiles  of  the 
uses  of  this  integrator.  Our  aim  is  two-fold;  to  interest  mathematicians 
so  that  this  method  may  be  placed  on  a firmer  mathematical  foundation 
and  to  inform  potential  users  so  that  they  may  appls'  this  algorithm  to 
their  particular  problem. 


TWO  BASIC  PUOBl.T.MS 

Before  we  discuss  the  sketch  and  the  examples  let  us  briefly  review 
two  problems  that  are  b.isic  to  any  ninnerical  solution  to  OI)i;S,  stiff  or 
not.  These  ]iroblems  are  truncation  error  and  stability.  According  to 
I'ahlquist  and  B jbrck  , ^ " (t  runcat  ion  errors')  are  committed  when  a 
liii)iti)ig  process  is  broken  off  before  one  has  come  to  the  limit  i))g 
value."  I'riuicat  ion  errors  result  from  mathematical  apjiroximat  ions . lor 
exanqile  they  arise  when  a finite  series  approximates  an  infinite  series 
or  when  a linear  function  iipjiroximates  a non-linear  one. 

Stability,  or  its  better  known  o(iposite,  instabi  1 i t>' . is  .issociateil 
with  the  idea  of  feedback.^  As  the  name  iinjilies,  ]iart  of  a program  or 
code  has  a loo[i  in  which  the  numbers  [iroduced  at  the  output  of  one  cycle 
are  used  as  the  input  for  the  next  cycle.  The  errors  associated  with  these 
numbers  may  then  be  amiilified  in  such  a way  as  to  destroy  the  solution. 

The  [)urpose  in  recalling  these  nemeses  is  to  Justify  the  effort 
that  has  gone  into  the  K-method  algorithm  to  reduce  truncation  erroi-  and 
to  maximize  stability  consistent  with  a reasonable  execution  time. 


^Numerical  Methods,  by  ti.  I):ihl((uist  and  A.  Bjbrck,  Trans.  b> 
Anderson,  11)74,  Trent  icc'-llall,  Inc.,  T.nglewood  Cliffs,  N.l,  p.  dj, 

“.See  for  example  Nuiiier i ca  1 Methods  for  Scientists  and  engineers,  by 
U.  W.  Hamming,  2nd  Iklition,  l‘)7.S,  McCrawuHill,  |nc.  p.  .1. 
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Higure  1 schematically  shows  the  main  functional  steps  in  the  K- 
method  algorithm.^  A third-order  predictor-corrector  method  is  emploj'ed. 
As  soon  as  the  initial  corrector  is  formed,  the  diagonal  of  the  dacohian 
(i.e.  3y.'/9y)  is  exaaiined  to  determine  those  dependent  variahles  that 
are  stif^  and*'those  that  are  not  (not  shown  in  Hig.  1).  This  is  done  in 
order  to  select  the  method  with  the  least  computational  overhead  for 
updating  each  of  the  initial  predictor  values. 

At  this  stage  of  the  algorithm  neither  the  predictor  values  nor 
corrector  values  are  presumed  acceptable  and  an  error  vector  is  gener- 
ated from  their  difference.  This  vector,  in  conjunction  with  the  Jacobian 
including,  now,  the  off-diagonal  elements,  is  used  to  modify  the  pre- 
dictor values.  The  modi f ications  to  the  predicted  values,  wliich  involve 
a matrix  inversion,  are  first  attempted  iteratively  using  a Gauss-Seidel 
method.  During  the  iteration,  checks  are  made  on  the  estimated  cora]iuta- 
tional  overhead  burden.  Should  the  iteration  method  prove  too  tedious, 
the  remaining  "non-converged"  correction  elements  are  solved  by  direct 
matrix  inversion,  llie  corrector  is  recomputed  and  another  error  vector 
is  generated. 

Truncation  error  is  then  checked  by  comparing  the  fourth  derivative 
of  each  dependent  variable  against  a predetermined  relative  error  toler- 
ance. If  any  one  of  these  variables  fails  tliis  test,  the  truncation 
error  is  judged  "poor,"  tlte  step  size,  h,  is  reduced  by  a factor  of  two 
and  this  cycle  of  the  computation  is  begun  anew.  This  test  is  especially 
useful,  as  we  shall  see  later,  in  the  case  of  discontinuous  driving 
functions. 

Wlien  all  the  dependent  variables  have  passed  the  truncation  test, 
(i.e.,  are  "OK"  in  Hig.  1)  a check  is  made  on  the  predictor-corrector 
agreement.  Should  all  elements  of  the  predictor-corrector  difference 
vector  be  less  than  a predetermined  minimum  error  (i.e.  are  "0K"1 , tlu 
corrector  values  are  accepted  as  the  solutions  at  this  time  step  and  tlie 
step  size  is  adjusted  for  the  next  cycle.  If  any  element  of  this  differ- 
ence vector  is  greater  than  a predetermined  maximum  error  tolerance  the 
agreement  is  judged  "poor,"  the  step  size  is  reduced  by  a factor  of  two 
.jnd  this  cycle  begun  anew.  Hor  the  intermediate  case,  in  which  all 
elements  of  this  difference  vector  are  less  than  the  maximum  error  test, 
and  in  which  at  least  one  is  greater  than  the  minimum  error  test,  the 
difference  vector  is  judged  "so-so."  When  this  is  the  case  the  ju-edicted 
values  are  again  modified  and  the  process  begun  anew  until  an  acceptable 
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An  earlier  version  has  been  reported.  Sec  "Description  and  Comparison 
of  the  K-Method  for  Performing  Numerical  Integration  of  Stiff  Onlinary 
Differential  H.quations,"  by  M.  1).  Kregel  and  1..  I..  l.ortie,  hill  Report 
No.  17.33,  July  1974,  AD  #A()03855. 
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Figure  1.  Schematic  of  K-method  of  integration;  h is  the  step  size. 
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solution  is  obtained.  It  is  important  to  note  that  the  K-method  conserves 
both  charge,  if  any,  and  chemical  balance  to  within  the  round-off  error 
of  the  machine  since  the  corrector  formulation  itself  is  conservative. 

The  methodology  outlined  in  Figure  1 and  above  shows  that  a solution 
over  one  time  step  is  considered  valid  if  the  corrected  values  map  into 
the  predicted  values  within  a specified  error  tolerance.  Since  small 
differences  in  the  predicted  values  are  magnified  by  a factor  of  the 
order  of  the  stiffness,  we  have  attempted  to  formulate  a predictor- 
corrector  scheme  which  has  minimum  error  and  maximum  stability.  How 
this  has  been  empirically  achieved  is  outlined  below. 

To  be  solved  are  vector  equations  which  may  be  cast  in  the  form  of 

Y'  = F - RY,  (1) 

where  Y is  the  dependent  vector,  and  a function  of  the  independent  variable 
-X,  F the  formation  vector  [ = F(Y,X)]  and  R the  removal  vector  [ = R(Y,XjJ. 
The  predictor  is  chosen  to  be  quadratic  in  form;  i.e., 

Y^^  = Yq  + A(X^  - Xp)  + B(Xj  - Xp)^  (2) 

where  (X  - X^)  is  the  local  step  size  and  where  the  subscript  zero 
denotes  the  current  location.  A and  B are  to  be  determined.  To  do  this 
we  require  two  independent  equations.  One  is  obtained  from  equation  121 

by  "looking  backward"  from  the  current  location,  X„,  to  the  time  X . 

P . u -1 

In  this  case  Y^  is  replaced  by  Y ^ which  has  been  evaluated,  i.e., 

Y-i  = Vq  + A(X_^  - Xq)  + B(X_^  - Xp)^  . (3) 


Obviously  another  equation  could  be  written  for  Y ^ but  an  alternate,  more 
stable  and  error-free  method  has  been  found.  Consider  the  derivative  of 
equation  (2)  at  Xj^ , namely, 

Yj*”  = A + 2B(Xj  - X(^)  . (4) 

and  the  predictor  in  the  form  of  equation  (1), 


P P 

Rl  Y^ 


(S1 


p 

Substituting  the  definition  of  Y^  from  equation  (2)  into  equation  (.S) 
and  equating  the  right  hand  sides  of  equation  (4)  and  equation  (5)  wc 
have 


A ^ 2B(X^  . x^)  = - R^PfYp  ^ A(X^  - X,,)  ^ B(Xj  - X^)"l. 


P P 

Provided  and  R^  are  known,  equation  (6)  provides  the  second 
required  to  determine  A and  B. 


((>) 

equat i on 
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1) 

l.et  us  now  consider  how  1'^  and  Rj  are  determined.  l'ij;ure  2 shows 
the  current  and  three  previous  discrete  values  of  the  formation  el^nent 
for  a ijiven  dependent  variable.  A parabola  of  the  fonii 

IHX]  = + C,(X  - - Xj,r’  , (7) 

is  passed  through  these  four  points.  llie  C.'s  (j  = l,2,.'5)  of  equation  (7) 
are  found  from  a least  squares  fit  where  th^  function  to  he  minimized 
with  respect  to  the  C.  is 


I W.  [f.  - F(X.)|^  - .1  ^ 1 < 0 . 


The  W.  are  weighting  functions,  chosen  so  that  periodic  fluctuations  in 

^ P 

the  1'^  values  will  not  propagate  into  . (The  relative  weights  are 

determined  by  adding  a quantity  a(-l)*  to  each  of  the  four  I . and  requir- 

P P * 

ing  that  the  identical  F.  be  found.)  R.  is  found  in  a similar  fashion. 
P P ^ ^ 

Once  F and  R,  are  determined,  A and  B can  be  tletermined,  and  conse- 

^ P ^ 
quently,  . 

The  corrector  is  given  by 

wliere  1)  and  B are  preselected  to  minimize  both  truncation  error  and 
noise  amplification,  and  to  maximize  relative  stability.  Values  for  the 
coefficients  in  equation  (9)  can  be  found  in  Reference  .1. 


liXAMPl.HS 

We  shall  now  consider  three  examples  of  the  types  of  problems  the 
K-method  has  been  called  upon  to  handle.  'I'hey  are  all  drawn  from  the 
field  of  acronomy.  Figure  3 shows  a linear  plot  of  a i)iecewisi‘  con- 
tinuous driving  function  [Q(t),  solid  line).  The  discontinuities  ai'e 
instantaneous  since  they  were  formed  by  reading  from  a HAl'A  block  with 
(.1)  and  (,J  + 1)  subscripts  interchanged.  (The  reverse  of  each  of  the 
slopes  in  Figure  7>  gives  the  desired  driving  function,  which  is  still 
discontinuous  in  the  first  derivative.)  The  dashed  linos  are  the 
resjionse  of  the  electron  density  and  a primary  positive  ion  densitv  , 
here  the  nitric  oxide  ion,  as  a function  of  time.  The  curves  have  been 
vertically  displaced  for  ease  in  reading.  It  is  seen  that  the  k-method 
enables  the  dependent  variables  to  follow  the  iitput  discontinuities  of 
the  driving  term.  (Departures  from  a perfect  m.itching  of  the  input 
slopes  can  be  explained  by  chemistry  competing  with  the  driving  term.) 
This  exiunple  demonstrates  the  effect  of  careful  monitoring  of  the  trun- 
cation error. 


Figure  2.  Schematic  of  a parabolic  fit  to  the  values  of  F at  the  current  and  three 

orevious  times  for  one  of  the  dependent  variables.  is  found  by  continuing 

this  parabola  to  the  time  X,. 


The  second  example,  Figure  4,  shows  histograms  of  the  number  of 
species  that  lie  in  the  decade  interval,  h/r.,  where  h is  the  local  steji 
size  and  where  t.  is  the  instantaneous  characteristic  time  constant  of 
the  ith  species.  (The  total  area  under  each  histogram  corresponds  to  64 
species.)  Tlie  numbers  to  the  far  right  are  the  decade  model  times  in 
seconds  (i.e.,  computer  execution  time  runs  from  bottom  to  top  in  this 
figure) . The  histograms  are  divided  into  a stiff  segment  to  the  right 
of  tlie  da|hed  line,  and  a non-stiff  segment  to  the  left.  On  the  first 
plot  (10  seconds)  only  a fe^  species  are  stiff  while  at  the  upper 
limit  of  this  integration  (10  seconds)  the  number  has  significantly 
increased,  with  stiffness  factors  (h/r^)  greater  than  10^. 

The  last  example  is  shown  in  Figures  5 and  6.  Figure  5 shows  the 
log  of  the  input  or  driving  function,  q(e),  plotted  against  the  log  of 

-4  -2 

time.  In  the  interval  10  - 10  seconds  those  negatively  charged 

species  more  strongly  coupled  to  the  driving  function,  (e.g.  e,  0.,  or 

Oj  ) follow  the  discontinuities  exhibited  by  the  driving  term.  These 

details  tend  to  be  "washed-out"  in  the  species  that  are  weakly  coupled 

2 

to  the  driving  term  (e.g.  CO^  or  CO^  ).  Near  10  seconds  the  strongly 

coupled  species  again  track  the  discontinuities  in  q(ej . The  d>Tiamic 

range  in  the  dependent  variables  is  about  six  orders  of  magnitude. 

Figure  (6)  shows  the  broad  range  response  of  the  neutral  species.  I'liis 

1 

graph  shows  those  species  that  follow  the  driving  term  [e.g.  0(  D)  , Nl-l))], 
those  that  are  independent  of  the  driving  term  [e.g.  N2O,  CO)  and  those 
that  tend  to  be  chemistry  dominated  [e.g.  H., , UNO,,  N.,0j|. 

In  summary,  we  have  empirically  derived  a third  order,  variable 
step  size  method  for  efficiently  handling  stiff  ODFS  that  appears  to 
work  for  discontinuous  driving  functions.  We  anticipate  with  some 
further  work  that  this  method  will  be  put  on  a firmer  mathematical  foun- 
dat ion. 


1 I 


A 


Histograms  that  characterize  the  stiffness  of  chemical  species  found  in 
daytime  deionization  model  of  the  atmosphere  at  an  altitude  of  30  km. 
(See  text  for  details). 
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Figure  5.  Modeled  response  of  the  negative  ion  densities  to  the 
driving  function  q(e) . 
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